
clc;
clear all;
close all;

%Read radius data file
load('data.mat');
R=radius;
D=R*2;

maxD=max(R)*2;
minD=min(R)*2;
%NumSphere=size(R,2);

%Calculate D-value
binInterval=0.3;
tempBin=[-10:binInterval:10];
maxDLog10=log10(maxD);
maxDLog10=min(find(tempBin>maxDLog10));
minDLog10=log10(minD);
minDLog10=max(find(tempBin<minDLog10));

phiBin=seqreverse([tempBin(minDLog10):binInterval:tempBin(maxDLog10)]);
numPhiBin=size(phiBin,2)-1;
Dref=1;
DBin=Dref*power(10,phiBin);

cumulCount=0;

for i=1:numPhiBin
    iBin=find(D<=DBin(i) & D>DBin(i+1));

        count=size(iBin,2);
    cumulCount=cumulCount+count;
    cumulFreq(i)=cumulCount;
    if count==0
        diameter(i)=(DBin(i)+DBin(i+1))/2;
    else
        diameter(i)=median(D(iBin));
    end;

    tempCenter=center(iBin,:);
    tempArea=area(iBin);
    bigR(i,:)=CalcBigR(tempCenter,totalArea,tempArea);
    %bigR(i,4)=diameter(i);
end;

%plot(bigR(:,2),bigR(:,1),'.-');
logDiameter=log10(diameter);
plot(logDiameter,bigR(:,1),'.-');

